#ifndef PCPS_COEFF_VALUE_HEADER
#define PCPS_COEFF_VALUE_HEADER

#include <src/pcps.h>
#include <src/input.h>
#include <src/correlator.h>
#include <src/hamiltonian.h>
#include <src/wavefunction.h>
#include <src/arena.h>
#include <src/blas_lapack.h>
#include <src/agp.h>
#include <src/pfaffian.h>
#include <src/mantissa_exponent.h>

namespace pcps {

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   computes the combined weights of the given correlators for the given configuration
  ///
  /// \param[in]      configs_per_site   The number of configurations on a single lattice site.
  /// \param[in]      correlators        The correlators.
  /// \param[in]      config             The lattice configuration.
  ///
  /// \return the combined weights of the given correlators for the given configuration
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template <class SCALAR>
  inline SCALAR correlator_weight(const int configs_per_site,
                                  const std::vector< pcps::Correlator<SCALAR> > & correlators,
                                  const int * const config) {

    // initialize the weight
    SCALAR retval = pcps::unity<SCALAR>();

    // compute each correlator's contribution
    for (int i = 0; i < correlators.size(); i++)
      retval *= correlators[i][config];

    // return the weight
    return retval;

  }

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   computes the combined weights of the given correlators for the given configuration
  ///
  /// \param[in]      configs_per_site   The number of configurations on a single lattice site.
  /// \param[in]      correlators        The correlators.
  /// \param[in]      config             The lattice configuration.
  /// \param[in,out]  weight             The wavefunction weight to which we add the correlators' contributions.
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template <class SCALAR>
  inline void correlator_weight(const int configs_per_site,
                                const std::vector< pcps::Correlator<SCALAR> > & correlators,
                                const int * const config,
                                pcps::MantExp<SCALAR> & weight) {

    // multiply the weight by each correlator's contribution
    for (int i = 0; i < correlators.size(); i++) {
      weight *= correlators[i][config];
    }

  }

  ///////////////////////////////////////////////////////////////////////////////////////////////////
  /// \brief   computes the combined weights of the given correlators for the given configuration
  ///
  /// \param[in]      configs_per_site   The number of configurations on a single lattice site.
  /// \param[in]      correlators        The correlators.
  /// \param[in]      config             The lattice configuration.
  ///
  /// \return the combined weights of the given correlators for the given configuration
  ///
  ///////////////////////////////////////////////////////////////////////////////////////////////////
  template <class SCALAR>
  inline SCALAR correlator_weight(const int configs_per_site,
                                  const std::vector< pcps::Correlator<SCALAR> *> & correlators,
                                  const int * const config) {

    // initialize the weight
    SCALAR retval = pcps::unity<SCALAR>();

    // compute each correlator's contribution
    for (int i = 0; i < correlators.size(); i++)
      retval *= (*correlators[i])[config];

    // return the weight
    return retval;

  }

  //-------------------------------------------------------------------------------
  // pcps::CoeffValueFunc -- A function object that computes the coefficient of a
  //                         given configuration.
  //
  //                         This object is templated on the scalar (SCALAR),
  //                         particle (PART), and reference (REF) types.
  //
  //                         The default of the template simply raises an error
  //                         to tell the user their choice of SCALAR, PART, and REF
  //                         was not valid.
  //
  //                         The various specializations of the template define
  //                         the actual behavior.
  //-------------------------------------------------------------------------------

  template <class SCALAR, class PART, class REF> class CoeffValueFunc {

    public:

    // constructor
    CoeffValueFunc(const pcps::Input & userinp,
                   const pcps::Hamiltonian<SCALAR> & ham,
                   const pcps::Wavefunction<SCALAR> & wfn,
                   pcps::Arena & arena) {
      CoeffValueFunc<SCALAR, PART, REF>::throw_error();
    }

    // function to initialize the internal data for the given configuration
    void init(const int * const config) {
      CoeffValueFunc<SCALAR, PART, REF>::throw_error();
    }

    // function for copying the state from another CoeffValueFunc
    void copy(const pcps::CoeffValueFunc<SCALAR, PART, REF> & other) {
      CoeffValueFunc<SCALAR, PART, REF>::throw_error();
    }

    // function to get the coefficient ratio for two configurations (new_config in numberator), which differ at the sites in diff_sites
    SCALAR get_coeff_ratio(const int n_diff_sites,
                           const int * const diff_sites,
                           const int * const new_config,
                           const int * const old_config) const {
      CoeffValueFunc<SCALAR, PART, REF>::throw_error();
      return pcps::zero<SCALAR>();
    }

    // function to update the internal data for a new configuration, which differs from the old configuration at the sites in diff_sites
    void update(const int n_diff_sites,
                const int * const diff_sites,
                const int * const new_config,
                const int * const old_config) {
      CoeffValueFunc<SCALAR, PART, REF>::throw_error();
    }

    // Function to get the ratio of the current coefficient's derivatives to the coefficient's value.
    // Derivatives with respect to each wavefunction variable are considered, although derivatives
    // with respect to most correlator elements are zero, so only non-trivial derivatives are returned.
    void get_compact_derivatives(const int * const config,
                                 int    * const cor_der_indices,
                                 SCALAR * const cor_der,
                                 SCALAR * const agp_der,
                                 SCALAR * const orb_der,
                                 SCALAR * const par_der) const {
      CoeffValueFunc<SCALAR, PART, REF>::throw_error();
    }

    // function to return the value of the wavefunction coefficient for the configuration that the object was last initialized for
    SCALAR initial_value() const {
      CoeffValueFunc<SCALAR, PART, REF>::throw_error();
      return pcps::zero<SCALAR>();
    }

    ///////////////////////////////////////////////////////////////////////////////////////////////////
    /// \brief   computes the wavefunction value for the given configuration
    ///
    /// \param[in]      config     The lattice configuration whose wavefunction value is desired.
    ///
    /// \return the value of the wavefunction in mantissa/exponent format.
    ///
    ///////////////////////////////////////////////////////////////////////////////////////////////////
    pcps::MantExp<SCALAR> get_wfn_value(const int * const config) {
      CoeffValueFunc<SCALAR, PART, REF>::throw_error();
      return pcps::MantExp<SCALAR>();
    }

    // function to tell the user that their choice of template parameters was invalid
    static void throw_error() {
      throw pcps::Exception("no implementation for pcps::CoeffValueFunc with template parameters SCALAR = %s, PART = %s, REF = %s")
              % pcps::type_lookup(typeid(SCALAR)) % pcps::type_lookup(typeid(PART)) % pcps::type_lookup(typeid(REF));
    }

  };

  //-------------------------------------------------------------------------------
  // pcps::CoeffValueFunc -- Specialization for the uniform reference.
  //-------------------------------------------------------------------------------

  template <class SCALAR, class PART> class CoeffValueFunc<SCALAR, PART, pcps::Ref_Uniform> {

    // private data
    private:
      const pcps::Wavefunction<SCALAR> & _wfn;
      const int _configs_per_site;
      const bool _get_cor_der;
      bool _initialized;
      SCALAR _init_val;

    public:

    // constructor
    CoeffValueFunc(const pcps::Input & userinp,
                   const pcps::Hamiltonian<SCALAR> & ham,
                   const pcps::Wavefunction<SCALAR> & wfn,
                   pcps::Arena & arena)
      : _wfn(wfn),
        _configs_per_site(userinp.configs_per_site()),
        _get_cor_der(userinp.optimize_corrs()),
        _initialized(false)
    {}
                   
    // function to initialize the internal data for the given configuration
    void init(const int * const config) {
      _init_val = pcps::correlator_weight(_configs_per_site, _wfn.correlators(), config);
      _initialized = true;
    }

    ///////////////////////////////////////////////////////////////////////////////////////////////////
    /// \brief   computes the wavefunction value for the given configuration
    ///
    /// \param[in]      config     The lattice configuration whose wavefunction value is desired.
    ///
    /// \return the value of the wavefunction in mantissa/exponent format.
    ///
    ///////////////////////////////////////////////////////////////////////////////////////////////////
    pcps::MantExp<SCALAR> get_wfn_value(const int * const config) {
      pcps::MantExp<SCALAR> retval(pcps::unity<SCALAR>());
      pcps::correlator_weight(_configs_per_site, _wfn.correlators(), config, retval);
      return retval;
    }

    // function to return the value of the wavefunction coefficient for the configuration that the object was last initialized for
    SCALAR initial_value() const { return _init_val; }

    // function for copying the state from another CoeffValueFunc
    void copy(const pcps::CoeffValueFunc<SCALAR, PART, pcps::Ref_Uniform> & other) {
      if (this == &other) return;
      _initialized = other._initialized;
    }

    // function to get the coefficient ratio for two configurations (new_config in numberator), which differ at the sites in diff_sites
    SCALAR get_coeff_ratio(const int n_diff_sites,
                           const int * const diff_sites,
                           const int * const new_config,
                           const int * const old_config) const {

      assert(_initialized);
      assert( n_diff_sites == 0 || n_diff_sites == 1 || n_diff_sites == 2 || n_diff_sites == 4 );

      // initialize the coefficient
      SCALAR retval = pcps::unity<SCALAR>();

      if (n_diff_sites == 0) {

      } else if (n_diff_sites == 1) {

        assert( typeid(PART) == typeid(pcps::Part_Spin) );

        // compute the correlator contribution
        retval *= pcps::correlator_weight(_configs_per_site, _wfn.touching_correlators(diff_sites[0]), new_config);
        retval /= pcps::correlator_weight(_configs_per_site, _wfn.touching_correlators(diff_sites[0]), old_config);

      } else if (n_diff_sites == 2) {

        // compute the correlator contribution
        retval *= pcps::correlator_weight(_configs_per_site, _wfn.touching_correlators(diff_sites[0], diff_sites[1]), new_config);
        retval /= pcps::correlator_weight(_configs_per_site, _wfn.touching_correlators(diff_sites[0], diff_sites[1]), old_config);

        // compute the sign from the assumed operation of a creation and destruction operator
        if (typeid(PART) == typeid(pcps::Part_Fermion) || typeid(PART) == typeid(pcps::Part_SpinlessFermion)) {
          if (typeid(PART) == typeid(pcps::Part_Fermion))
            assert(diff_sites[0] % 2 == diff_sites[1] % 2);
          const int I = std::min(diff_sites[0], diff_sites[1]);
          const int J = std::max(diff_sites[0], diff_sites[1]);
          for (int i = I+1; i < J; i++) {
            assert( old_config[i] == 0 || old_config[i] == 1 );
            retval *= 1 - 2 * old_config[i];
          }
        }

      } else if (n_diff_sites == 4) {

        if (typeid(PART) != typeid(pcps::Part_Spin))
          throw pcps::Exception( (boost::format("error in get_coeff_ratio for uniform reference.  n_diff_sites = 4 is allowed only for particles of type \"%s\"")
                                  % pcps::keyword::Part_Spin).str() );

        // compute the correlator contribution
        const int * const ds = diff_sites;
        retval *= pcps::correlator_weight(_configs_per_site, _wfn.touching_correlators(ds[0], ds[1], ds[2], ds[3]), new_config);
        retval /= pcps::correlator_weight(_configs_per_site, _wfn.touching_correlators(ds[0], ds[1], ds[2], ds[3]), old_config);

      } else {

        throw pcps::Exception("error in get_coeff_ratio for uniform reference.  n_diff_sites must be 0, 1, 2, or 4.");

      }

      return retval;

    }

    // function to update the internal data for a new configuration, which differs from the old configuration at the sites in diff_sites
    void update(const int n_diff_sites,
                const int * const diff_sites,
                const int * const new_config,
                const int * const old_config) {
      assert(_initialized);
    }

    // Function to get the ratio of the current coefficient's derivatives to the coefficient's value.
    // Derivatives with respect to each wavefunction variable are considered, although derivatives
    // with respect to most correlator elements are zero, so only non-trivial derivatives are returned.
    void get_compact_derivatives(const int * const config,
                                 int    * const cor_der_indices,
                                 SCALAR * const cor_der,
                                 SCALAR * const agp_der,
                                 SCALAR * const orb_der,
                                 SCALAR * const par_der) const {

      // compute correlator derivatives (skip fixed correlators)
      if (_get_cor_der) {
        int c = 0;
        for (typename std::vector<pcps::Correlator<SCALAR> >::const_iterator corr = _wfn.correlators().begin(); corr != _wfn.correlators().end(); corr++)
          if (!corr->fixed()) {
            int cmpd = 0;
            const int * sites = &corr->site_vec()[0];
            const int nsites = corr->site_vec().size();
            for (int i = 0; i < nsites; i++)
              cmpd = cmpd * _configs_per_site + config[sites[i]];
            cor_der_indices[c] = cmpd;
            cor_der[c] = pcps::unity<SCALAR>() / (*corr)[cmpd];
            c++;
          }
      }

      // there are no agp weight, orbital coefficient, or pairing matrix derivatives for the uniform reference

    }

  };

  //-------------------------------------------------------------------------------
  // pcps::CoeffValueFunc -- Specialization for fermions and the restricted
  //                         antisymmetrized geminal power reference.
  //-------------------------------------------------------------------------------

  template <class SCALAR> class CoeffValueFunc<SCALAR, pcps::Part_Fermion, pcps::Ref_AGP> {

    // protected data
    protected:
      const pcps::Wavefunction<SCALAR> & _wfn;  ///< reference to the wavefunction object
      const int _configs_per_site;              ///< size of one site's hilbert space
      const int _nap;                           ///< number of alpha particles (assumed to be equal to number of beta)
      const int _nas;                           ///< number of alpha sites (assumed to be equal to number of beta)
      const bool _get_cor_der;                  ///< whether to compute derivative ratios w.r.t. correlator variables
      const bool _get_agp_der;                  ///< whether to compute derivative ratios w.r.t. agp weight variables
      const bool _get_orb_der;                  ///< whether to compute derivative ratios w.r.t. restricted orbital coefficient variables
      const bool _get_par_der;                  ///< whether to compute derivative ratios w.r.t. the pairing matrix elements
      const bool _use_pairing_mat;              ///< whether to build the pairing matrix gamma from the wavefunction's pairing matrix or agp variables
      bool _initialized;                        ///< whether or not the object has been initialized with a lattice configuration
      SCALAR _init_val;                         ///< the value of the wavefunction for the most recently initialized configuration
      SCALAR * _gamma;                          ///< matrix of the possible elements of the agp matrix
      SCALAR * _agp_mat;                        ///< workspace used for constructing and inverting the agp matrix
      SCALAR * _inv_mat;                        ///< the inverse of the agp matrix
      SCALAR * _wv1;                            ///< workspace
      SCALAR * _wv2;                            ///< workspace
      SCALAR * _wv3;                            ///< workspace
      int    * _a_sites;                        ///< the currently occupied alpha sites
      int    * _b_sites;                        ///< the currently occupied beta sites
      int    * _a_sites_w;                      ///< workspace for alpha sites
      int    * _b_sites_w;                      ///< workspace for beta sites
      int    * _ipiv;                           ///< workspace for recording the pivoting that occurs in the LU decomposition
      int    * _inv_perm;                       ///< workspace for inverting the permutations associated with pivoting

    public:

    // constructor
    CoeffValueFunc(const pcps::Input & userinp,
                   const pcps::Hamiltonian<SCALAR> & ham,
                   const pcps::Wavefunction<SCALAR> & wfn,
                   pcps::Arena & arena,
                   const bool use_pairing_mat = false)
      : _wfn(wfn),
        _configs_per_site(userinp.configs_per_site()),
        _nap(userinp.nparticles()/2),
        _nas(userinp.nsites()/2),
        _get_cor_der(userinp.optimize_corrs()),
        _get_agp_der(userinp.optimize_agp_weights()),
        _get_orb_der(userinp.optimize_orb_coeffs()),
        _get_par_der(userinp.optimize_pairing_matrix()),
        _use_pairing_mat(use_pairing_mat),
        _initialized(false)
    {
      const SCALAR z = pcps::zero<SCALAR>();
      _gamma     = arena.allocate_array<SCALAR>( _nas * _nas, z );
      _agp_mat   = arena.allocate_array<SCALAR>( _nap * _nap, z );
      _inv_mat   = arena.allocate_array<SCALAR>( _nap * _nap, z );
      _wv1       = arena.allocate_array<SCALAR>( _nap       , z );
      _wv2       = arena.allocate_array<SCALAR>( _nap       , z );
      _wv3       = arena.allocate_array<SCALAR>( _nap       , z );
      _a_sites   = arena.allocate_array<int   >( _nap       , 0 );
      _b_sites   = arena.allocate_array<int   >( _nap       , 0 );
      _a_sites_w = arena.allocate_array<int   >( _nap       , 0 );
      _b_sites_w = arena.allocate_array<int   >( _nap       , 0 );
      _ipiv      = arena.allocate_array<int   >( _nap       , 0 );
      _inv_perm  = arena.allocate_array<int   >( _nap       , 0 );
    }

    // function to initialize the internal data for the given configuration
    void init(const int * const config) {

      // initialize the agp inverse matrix, gamma matrix, and alpha and beta site lists and get the initial determinant value
      pcps::agp::init_matrices(_wfn, _nap, _nas, config, _init_val, _a_sites, _b_sites, _gamma, _inv_mat, _agp_mat, _ipiv, _inv_perm, _use_pairing_mat);

      // get correlator contribution to initial value
      _init_val *= pcps::correlator_weight(_configs_per_site, _wfn.correlators(), config);

      // set the initialized flag to true
      _initialized = true;

    }

    ///////////////////////////////////////////////////////////////////////////////////////////////////
    /// \brief   computes the wavefunction value for the given configuration
    ///
    /// \param[in]      config     The lattice configuration whose wavefunction value is desired.
    ///
    /// \return the value of the wavefunction in mantissa/exponent format.
    ///
    ///////////////////////////////////////////////////////////////////////////////////////////////////
    pcps::MantExp<SCALAR> get_wfn_value(const int * const config) {

      // initialize the wavefunction value
      pcps::MantExp<SCALAR> retval(pcps::unity<SCALAR>());

      // get the correlators' contribution
      pcps::correlator_weight(_configs_per_site, _wfn.correlators(), config, retval);

      // get the pairing matrix's contribution
      retval *= pcps::agp::get_det_value(_wfn, _nap, _nas, config, _gamma, _a_sites_w, _b_sites_w, _agp_mat, _ipiv, _use_pairing_mat);

      // return the wavefunction value
      return retval;

    }

    // function to return the value of the wavefunction coefficient for the configuration that the object was last initialized for
    SCALAR initial_value() const { return _init_val; }

    // function for copying the state from another CoeffValueFunc
    void copy(const pcps::CoeffValueFunc<SCALAR, pcps::Part_Fermion, pcps::Ref_AGP> & other) {
      if (this == &other) return;
      _initialized = other._initialized;
      std::memcpy( _gamma    , other._gamma    , sizeof(SCALAR) * size_t(_nas * _nas) );
      std::memcpy( _agp_mat  , other._agp_mat  , sizeof(SCALAR) * size_t(_nap * _nap) );
      std::memcpy( _inv_mat  , other._inv_mat  , sizeof(SCALAR) * size_t(_nap * _nap) );
      std::memcpy( _wv1      , other._wv1      , sizeof(SCALAR) * size_t(_nap       ) );
      std::memcpy( _wv2      , other._wv2      , sizeof(SCALAR) * size_t(_nap       ) );
      std::memcpy( _wv3      , other._wv3      , sizeof(SCALAR) * size_t(_nap       ) );
      std::memcpy( _a_sites  , other._a_sites  , sizeof(int   ) * size_t(_nap       ) );
      std::memcpy( _b_sites  , other._b_sites  , sizeof(int   ) * size_t(_nap       ) );
      std::memcpy( _ipiv     , other._ipiv     , sizeof(int   ) * size_t(_nap       ) );
      std::memcpy( _inv_perm , other._inv_perm , sizeof(int   ) * size_t(_nap       ) );
    }

    // function to get the coefficient ratio for two configurations (new_config in numberator), which differ at the sites in diff_sites
    SCALAR get_coeff_ratio(const int n_diff_sites,
                           const int * const diff_sites,
                           const int * const new_config,
                           const int * const old_config) const {

      // check that this object is initialized
      assert(_initialized);

      // check that the number of sites affected by the move is supported
      assert( n_diff_sites == 0 || n_diff_sites == 2 || n_diff_sites == 4 );

      // return 1 if no sites were changed
      if (n_diff_sites == 0)
        return pcps::unity<SCALAR>();

      // initialize the coefficient
      SCALAR retval = pcps::unity<SCALAR>();

      // compute the correlators' contribution to the new coefficient
      if (n_diff_sites == 2) {
        assert( diff_sites[0] % 2 == diff_sites[1] % 2 );
        retval *= pcps::correlator_weight(_configs_per_site, _wfn.touching_correlators(diff_sites[0], diff_sites[1]), new_config);
        retval /= pcps::correlator_weight(_configs_per_site, _wfn.touching_correlators(diff_sites[0], diff_sites[1]), old_config);
      } else if (n_diff_sites == 4) {
        const int * const ds = diff_sites;
        retval *= pcps::correlator_weight(_configs_per_site, _wfn.correlators(), new_config);
        retval /= pcps::correlator_weight(_configs_per_site, _wfn.correlators(), old_config);
        //retval *= pcps::correlator_weight(_configs_per_site, _wfn.touching_correlators(ds[0], ds[1], ds[2], ds[3]), new_config);
        //retval /= pcps::correlator_weight(_configs_per_site, _wfn.touching_correlators(ds[0], ds[1], ds[2], ds[3]), old_config);
      } else throw pcps::Exception( (boost::format("%i is not a valid value for n_diff_sites in pcps::CoeffValueFunc::get_coeff_ratio") % n_diff_sites).str() );

      // multiply by the agp determinant ratio
      retval *= pcps::agp::det_ratio(_nap, _nas, n_diff_sites, diff_sites, new_config, old_config, _gamma, _inv_mat, _a_sites, _b_sites, _wv1, _wv2, _wv3);

      // return the coefficient
      return retval;

    }

    // function to update the internal data for a new configuration, which differs from the old configuration at the sites in diff_sites
    void update(const int n_diff_sites,
                const int * const diff_sites,
                const int * const new_config,
                const int * const old_config) {

      // check that this object is initialized
      assert(_initialized);

      // check that the number of sites affected by the move is supported
      assert(n_diff_sites == 0 || n_diff_sites == 2 || n_diff_sites == 4);

      // there is nothing to do if no sites were affected
      if (n_diff_sites == 0) return;

      // to avoid walking through nodes, re-initialize for double moves
      if (n_diff_sites == 4) {
        this->init(new_config);
        return;
      }

      // update the inverse agp matrix and the alpha and beta site lists
      pcps::agp::update(_nap, _nas, n_diff_sites, diff_sites, new_config, old_config, _gamma, _inv_mat, _a_sites, _b_sites, _wv1, _wv2, _wv3);

    }

    // Function to get the ratio of the current coefficient's derivatives to the coefficient's value.
    // Derivatives with respect to each wavefunction variable are considered, although derivatives
    // with respect to most correlator elements are zero, so only non-trivial derivatives are returned.
    void get_compact_derivatives(const int * const config,
                                 int    * const cor_der_indices,
                                 SCALAR * const cor_der,
                                 SCALAR * const agp_der,
                                 SCALAR * const orb_der,
                                 SCALAR * const par_der) const {

      // compute correlator derivative ratios (skip fixed correlators)
      if (_get_cor_der) {
        int c = 0;
        for (typename std::vector<pcps::Correlator<SCALAR> >::const_iterator corr = _wfn.correlators().begin(); corr != _wfn.correlators().end(); corr++)
          if (!corr->fixed()) {
            int cmpd = 0;
            const int * sites = &corr->site_vec()[0];
            const int nsites = corr->site_vec().size();
            for (int i = 0; i < nsites; i++)
              cmpd = cmpd * _configs_per_site + config[sites[i]];
            cor_der_indices[c] = cmpd;
            cor_der[c] = pcps::unity<SCALAR>() / (*corr)[cmpd];
            c++;
          }
      }

      // compute derivative ratios for the agp weights and orbital coefficients
      pcps::agp::der_ratios(_wfn, _nap, _nas, _get_agp_der, _get_orb_der, _get_par_der, _gamma, _inv_mat, _a_sites, _b_sites, agp_der, orb_der, par_der);

    }

  };

  //-------------------------------------------------------------------------------
  // pcps::CoeffValueFunc -- Specialization for spins and the restricted
  //                         antisymmetrized geminal power reference.
  //-------------------------------------------------------------------------------

  template <class SCALAR> class CoeffValueFunc<SCALAR, pcps::Part_Spin, pcps::Ref_AGP> {

    // private data
    protected:
      const pcps::Wavefunction<SCALAR> & _wfn;  ///< reference to the wavefunction object
      const int _configs_per_site;              ///< size of one site's hilbert space
      const int _nap;                           ///< number of alpha particles (assumed to be equal to number of beta)
      const int _nas;                           ///< number of alpha sites (assumed to be equal to number of beta)
      const bool _get_cor_der;                  ///< whether to compute derivative ratios w.r.t. correlator variables
      const bool _get_agp_der;                  ///< whether to compute derivative ratios w.r.t. agp weight variables
      const bool _get_orb_der;                  ///< whether to compute derivative ratios w.r.t. restricted orbital coefficient variables
      const bool _get_par_der;                  ///< whether to compute derivative ratios w.r.t. the pairing matrix elements
      const bool _use_pairing_mat;              ///< whether to build the pairing matrix gamma from the wavefunction's pairing matrix or agp variables
      bool _initialized;                        ///< whether or not the object has been initialized with a lattice configuration
      SCALAR _init_val;                         ///< the value of the wavefunction for the most recently initialized configuration
      SCALAR * _gamma;                          ///< matrix of the possible elements of the agp matrix
      SCALAR * _agp_mat;                        ///< workspace used for constructing and inverting the agp matrix
      SCALAR * _inv_mat;                        ///< the inverse of the agp matrix
      SCALAR * _wv1;                            ///< workspace
      SCALAR * _wv2;                            ///< workspace
      SCALAR * _wv3;                            ///< workspace
      int    * _a_sites;                        ///< the currently occupied alpha sites
      int    * _b_sites;                        ///< the currently occupied beta sites
      int    * _a_sites_w;                      ///< workspace for alpha sites
      int    * _b_sites_w;                      ///< workspace for beta sites
      int    * _new_conf;                       ///< fictitious configuration of fermions representing the spins for use with agp
      int    * _old_conf;                       ///< fictitious configuration of fermions representing the spins for use with agp
      int    * _ferm_diff_sites;                ///< list of sites differing between the old and new fermion configurations
      int    * _ipiv;                           ///< workspace for recording the pivoting that occurs in the LU decomposition
      int    * _inv_perm;                       ///< workspace for inverting the permutations associated with pivoting

    public:

    // constructor
    CoeffValueFunc(const pcps::Input & userinp,
                   const pcps::Hamiltonian<SCALAR> & ham,
                   const pcps::Wavefunction<SCALAR> & wfn,
                   pcps::Arena & arena,
                   const bool use_pairing_mat = false)
      : _wfn(wfn),
        _configs_per_site(userinp.configs_per_site()),
        _nap(userinp.nsites()/2),
        _nas(userinp.nsites()),
        _get_cor_der(userinp.optimize_corrs()),
        _get_agp_der(userinp.optimize_agp_weights()),
        _get_orb_der(userinp.optimize_orb_coeffs()),
        _get_par_der(userinp.optimize_pairing_matrix()),
        _use_pairing_mat(use_pairing_mat),
        _initialized(false)
    {

      // check that we have spin-1/2 spins
      if ( _configs_per_site != 2 || std::abs( userinp.particle_spin() - 0.5 ) > 1.0e-6 )
        throw pcps::Exception( (boost::format("pcps::CoeffValueFunc for spins and agp reference requires spin 1/2 spins.")).str() );

      // check that we have an equal number of up and down spins
      if ( userinp.nsites() % 2 != 0 || std::abs( userinp.z_spin() ) > 1.0e-6 )
        throw pcps::Exception( (boost::format("pcps::CoeffValueFunc for spins and agp reference requires an equal number of up and down spins.")).str() );

      // allocate arrays
      const SCALAR z = pcps::zero<SCALAR>();
      _gamma           = arena.allocate_array<SCALAR>( _nas * _nas, z );
      _agp_mat         = arena.allocate_array<SCALAR>( _nap * _nap, z );
      _inv_mat         = arena.allocate_array<SCALAR>( _nap * _nap, z );
      _wv1             = arena.allocate_array<SCALAR>( _nap       , z );
      _wv2             = arena.allocate_array<SCALAR>( _nap       , z );
      _wv3             = arena.allocate_array<SCALAR>( _nap       , z );
      _a_sites         = arena.allocate_array<int   >( _nap       , 0 );
      _b_sites         = arena.allocate_array<int   >( _nap       , 0 );
      _a_sites_w       = arena.allocate_array<int   >( _nap       , 0 );
      _b_sites_w       = arena.allocate_array<int   >( _nap       , 0 );
      _new_conf        = arena.allocate_array<int   >( 2 * _nas   , 0 );
      _old_conf        = arena.allocate_array<int   >( 2 * _nas   , 0 );
      _ferm_diff_sites = arena.allocate_array<int   >( 4          , 0 );
      _ipiv            = arena.allocate_array<int   >( _nap       , 0 );
      _inv_perm        = arena.allocate_array<int   >( _nap       , 0 );

    }

    // function to convert the spin configuration to a fermion configuration
    void spin_conf_to_fermion_conf(const int * const sconfig, int * const fconfig) {

      // loop over the spin sites
      for (int i = 0; i < _nas; i++) {

        // zero the alpha and beta occupations of the ith site
        fconfig[2*i+0] = 0;
        fconfig[2*i+1] = 0;

        // set either alpha or beta as occupied
        fconfig[2*i+1-sconfig[i]] = 1;

      }

    }

    // function to initialize the internal data for the given configuration
    void init(const int * const config) {

      // convert the spin configuration to a fermion configuration
      this->spin_conf_to_fermion_conf(config, _old_conf);

      // initialize the agp inverse matrix, gamma matrix, and alpha and beta site lists and get the initial determinant value
      pcps::agp::init_matrices(_wfn, _nap, _nas, _old_conf, _init_val, _a_sites, _b_sites, _gamma, _inv_mat, _agp_mat, _ipiv, _inv_perm, _use_pairing_mat);

      // get correlator contribution to initial value
      _init_val *= pcps::correlator_weight(_configs_per_site, _wfn.correlators(), config);

      // set the initialized flag to true
      _initialized = true;

    }

    ///////////////////////////////////////////////////////////////////////////////////////////////////
    /// \brief   computes the wavefunction value for the given configuration
    ///
    /// \param[in]      config     The lattice configuration whose wavefunction value is desired.
    ///
    /// \return the value of the wavefunction in mantissa/exponent format.
    ///
    ///////////////////////////////////////////////////////////////////////////////////////////////////
    pcps::MantExp<SCALAR> get_wfn_value(const int * const config) {

      // initialize the wavefunction value
      pcps::MantExp<SCALAR> retval(pcps::unity<SCALAR>());

      // get the correlators' contribution
      pcps::correlator_weight(_configs_per_site, _wfn.correlators(), config, retval);

      // convert the spin configuration to a fermion configuration
      this->spin_conf_to_fermion_conf(config, _old_conf);

      // get the pairing matrix's contribution
      retval *= pcps::agp::get_det_value(_wfn, _nap, _nas, _old_conf, _gamma, _a_sites_w, _b_sites_w, _agp_mat, _ipiv, _use_pairing_mat);

      // return the wavefunction value
      return retval;

    }

    // function to return the value of the wavefunction coefficient for the configuration that the object was last initialized for
    SCALAR initial_value() const { return _init_val; }

    // function for copying the state from another CoeffValueFunc
    void copy(const pcps::CoeffValueFunc<SCALAR, pcps::Part_Spin, pcps::Ref_AGP> & other) {
      if (this == &other) return;
      _initialized = other._initialized;
      std::memcpy( _gamma          , other._gamma          , sizeof(SCALAR) * size_t(_nas * _nas) );
      std::memcpy( _agp_mat        , other._agp_mat        , sizeof(SCALAR) * size_t(_nap * _nap) );
      std::memcpy( _inv_mat        , other._inv_mat        , sizeof(SCALAR) * size_t(_nap * _nap) );
      std::memcpy( _wv1            , other._wv1            , sizeof(SCALAR) * size_t(_nap       ) );
      std::memcpy( _wv2            , other._wv2            , sizeof(SCALAR) * size_t(_nap       ) );
      std::memcpy( _wv3            , other._wv3            , sizeof(SCALAR) * size_t(_nap       ) );
      std::memcpy( _a_sites        , other._a_sites        , sizeof(int   ) * size_t(_nap       ) );
      std::memcpy( _b_sites        , other._b_sites        , sizeof(int   ) * size_t(_nap       ) );
      std::memcpy( _new_conf       , other._new_conf       , sizeof(int   ) * size_t(2 * _nas   ) );
      std::memcpy( _old_conf       , other._old_conf       , sizeof(int   ) * size_t(2 * _nas   ) );
      std::memcpy( _ferm_diff_sites, other._ferm_diff_sites, sizeof(int   ) * size_t(4          ) );
      std::memcpy( _ipiv           , other._ipiv           , sizeof(int   ) * size_t(_nap       ) );
      std::memcpy( _inv_perm       , other._inv_perm       , sizeof(int   ) * size_t(_nap       ) );
    }

    // function to get the coefficient ratio for two configurations (new_config in numberator), which differ at the sites in diff_sites
    SCALAR get_coeff_ratio(const int n_diff_sites,
                           const int * const diff_sites,
                           const int * const new_config,
                           const int * const old_config) {

      // check that this object is initialized
      assert(_initialized);

      // check that the number of sites affected by the move is supported
      assert( n_diff_sites == 0 || n_diff_sites == 2 );

      // return 1 if no sites were changed
      if (n_diff_sites == 0)
        return pcps::unity<SCALAR>();

      // convert the spin configurations to fermion configurations
      this->spin_conf_to_fermion_conf(new_config, _new_conf);
      this->spin_conf_to_fermion_conf(old_config, _old_conf);

      // initialize the coefficient
      SCALAR retval = pcps::unity<SCALAR>();

      if (n_diff_sites == 2) {

        // compute the correlators' contribution to the new coefficient
        retval *= pcps::correlator_weight(_configs_per_site, _wfn.touching_correlators(diff_sites[0], diff_sites[1]), new_config);
        retval /= pcps::correlator_weight(_configs_per_site, _wfn.touching_correlators(diff_sites[0], diff_sites[1]), old_config);

        // compile the list of fermion sites that differ (we assume the operator in question is RL)
        _ferm_diff_sites[0] = 2*diff_sites[0]+0;
        _ferm_diff_sites[1] = 2*diff_sites[0]+1;
        _ferm_diff_sites[2] = 2*diff_sites[1]+1;
        _ferm_diff_sites[3] = 2*diff_sites[1]+0;

        // multiply by the agp determinant ratio
        retval *= pcps::agp::det_ratio(_nap, _nas, 4, _ferm_diff_sites, _new_conf, _old_conf, _gamma, _inv_mat, _a_sites, _b_sites, _wv1, _wv2, _wv3);

      } else throw pcps::Exception( (boost::format("%i is not a valid value for n_diff_sites in pcps::CoeffValueFunc::get_coeff_ratio") % n_diff_sites).str() );

      // return the coefficient
      return retval;

    }

    // function to update the internal data for a new configuration, which differs from the old configuration at the sites in diff_sites
    void update(const int n_diff_sites,
                const int * const diff_sites,
                const int * const new_config,
                const int * const old_config) {

      // check that this object is initialized
      assert(_initialized);

      // check that the number of sites affected by the move is supported
      assert( n_diff_sites == 0 || n_diff_sites == 2 );

      // there is nothing to do if no sites were affected
      if (n_diff_sites == 0) return;

      // convert the spin configurations to fermion configurations
      this->spin_conf_to_fermion_conf(new_config, _new_conf);
      this->spin_conf_to_fermion_conf(old_config, _old_conf);

      // compile the list of fermion sites that differ (we assume the operator in question is RL)
      _ferm_diff_sites[0] = 2*diff_sites[0]+0;
      _ferm_diff_sites[1] = 2*diff_sites[0]+1;
      _ferm_diff_sites[2] = 2*diff_sites[1]+1;
      _ferm_diff_sites[3] = 2*diff_sites[1]+0;

      // update the inverse agp matrix and the alpha and beta site lists
      pcps::agp::update(_nap, _nas, 4, _ferm_diff_sites, _new_conf, _old_conf, _gamma, _inv_mat, _a_sites, _b_sites, _wv1, _wv2, _wv3);

    }

    // Function to get the ratio of the current coefficient's derivatives to the coefficient's value.
    // Derivatives with respect to each wavefunction variable are considered, although derivatives
    // with respect to most correlator elements are zero, so only non-trivial derivatives are returned.
    void get_compact_derivatives(const int * const config,
                                 int    * const cor_der_indices,
                                 SCALAR * const cor_der,
                                 SCALAR * const agp_der,
                                 SCALAR * const orb_der,
                                 SCALAR * const par_der) const {

      // compute correlator derivative ratios (skip fixed correlators)
      if (_get_cor_der) {
        int c = 0;
        for (typename std::vector<pcps::Correlator<SCALAR> >::const_iterator corr = _wfn.correlators().begin(); corr != _wfn.correlators().end(); corr++)
          if (!corr->fixed()) {
            int cmpd = 0;
            const int * sites = &corr->site_vec()[0];
            const int nsites = corr->site_vec().size();
            for (int i = 0; i < nsites; i++)
              cmpd = cmpd * _configs_per_site + config[sites[i]];
            cor_der_indices[c] = cmpd;
            cor_der[c] = pcps::unity<SCALAR>() / (*corr)[cmpd];
            c++;
          }
      }

      // compute derivative ratios for the agp weights and orbital coefficients
      pcps::agp::der_ratios(_wfn, _nap, _nas, _get_agp_der, _get_orb_der, _get_par_der, _gamma, _inv_mat, _a_sites, _b_sites, agp_der, orb_der, par_der);

    }

  };

  //-------------------------------------------------------------------------------
  // pcps::CoeffValueFunc -- Specialization for fermions and the pairing reference.
  //-------------------------------------------------------------------------------

  template <class SCALAR> class CoeffValueFunc<SCALAR, pcps::Part_Fermion, pcps::Ref_Pairing>
  : public CoeffValueFunc<SCALAR, pcps::Part_Fermion, pcps::Ref_AGP> {

    public:

    // constructor
    CoeffValueFunc(const pcps::Input & userinp,
                   const pcps::Hamiltonian<SCALAR> & ham,
                   const pcps::Wavefunction<SCALAR> & wfn,
                   pcps::Arena & arena)
      : CoeffValueFunc<SCALAR, pcps::Part_Fermion, pcps::Ref_AGP>(userinp, ham, wfn, arena, true)
    {}

  };

  //-------------------------------------------------------------------------------
  // pcps::CoeffValueFunc -- Specialization for spins and the pairing reference.
  //-------------------------------------------------------------------------------

  template <class SCALAR> class CoeffValueFunc<SCALAR, pcps::Part_Spin, pcps::Ref_Pairing>
  : public CoeffValueFunc<SCALAR, pcps::Part_Spin, pcps::Ref_AGP> {

    public:

    // constructor
    CoeffValueFunc(const pcps::Input & userinp,
                   const pcps::Hamiltonian<SCALAR> & ham,
                   const pcps::Wavefunction<SCALAR> & wfn,
                   pcps::Arena & arena)
      : CoeffValueFunc<SCALAR, pcps::Part_Spin, pcps::Ref_AGP>(userinp, ham, wfn, arena, true)
    {}

  };

  //-------------------------------------------------------------------------------
  // pcps::CoeffValueFunc -- Specialization for fermions and the pfaffian reference
  //-------------------------------------------------------------------------------

  template <class SCALAR> class CoeffValueFunc<SCALAR, pcps::Part_Fermion, pcps::Ref_Pfaffian> {

    // protected data
    protected:
      const pcps::Wavefunction<SCALAR> & _wfn;  ///< reference to the wavefunction object
      const int _configs_per_site;              ///< size of one site's hilbert space
      const int _np;                            ///< number of particles (assumed to be even)
      const int _ns;                            ///< number of sites
      const bool _get_cor_der;                  ///< whether to compute derivative ratios w.r.t. correlator variables
      const bool _get_agp_der;                  ///< whether to compute derivative ratios w.r.t. agp weight variables
      const bool _get_orb_der;                  ///< whether to compute derivative ratios w.r.t. restricted orbital coefficient variables
      const bool _get_par_der;                  ///< whether to compute derivative ratios w.r.t. the pairing matrix elements
      bool _initialized;                        ///< whether or not the object has been initialized with a lattice configuration
      int _update_count;                        ///< the number of times the object has been updated without being re-initialized
      pcps::MantExp<SCALAR> _init_val;          ///< the value of the wavefunction for the most recently initialized configuration
      SCALAR * _pm;                             ///< space to hold the full antisymmetric pairing matrix
      SCALAR * _inv_mat;                        ///< the inverse of the occupied pairing matrix
      SCALAR * _work_mat;                       ///< workspace
      SCALAR * _wv1;                            ///< workspace
      SCALAR * _wv2;                            ///< workspace
      SCALAR * _wv3;                            ///< workspace
      int    * _row_indices;                    ///< row indices of occupied sites in the occupied pairing matrix (-1 indicates not occupied)
      int    * _occ_sites;                      ///< the currently occupied sites (in the same order they appear in the occupied pairing matrix)
      int    * _ipiv;                           ///< workspace for recording the pivoting that occurs in the LU decomposition
      int    * _iwork;                          ///< integer workspace

    public:

    // constructor
    CoeffValueFunc(const pcps::Input & userinp,
                   const pcps::Hamiltonian<SCALAR> & ham,
                   const pcps::Wavefunction<SCALAR> & wfn,
                   pcps::Arena & arena)
      : _wfn(wfn),
        _configs_per_site(userinp.configs_per_site()),
        _np(userinp.nparticles()),
        _ns(userinp.nsites()),
        _get_cor_der(userinp.optimize_corrs()),
        _get_agp_der(userinp.optimize_agp_weights()),
        _get_orb_der(userinp.optimize_orb_coeffs()),
        _get_par_der(userinp.optimize_pairing_matrix()),
        _initialized(false),
        _update_count(0)
    {
      assert( _np % 2 == 0 );
      const SCALAR z = pcps::zero<SCALAR>();
      _pm          = arena.allocate_array<SCALAR>( _ns * _ns , z );
      _inv_mat     = arena.allocate_array<SCALAR>( _np * _np , z );
      _work_mat    = arena.allocate_array<SCALAR>( _np * _np , z );
      _wv1         = arena.allocate_array<SCALAR>( _np       , z );
      _wv2         = arena.allocate_array<SCALAR>( _np       , z );
      _wv3         = arena.allocate_array<SCALAR>( _np       , z );
      _row_indices = arena.allocate_array<int   >( _ns       , 0 );
      _occ_sites   = arena.allocate_array<int   >( _np       , 0 );
      _ipiv        = arena.allocate_array<int   >( _np       , 0 );
      _iwork       = arena.allocate_array<int   >( _ns       , 0 );
    }

    // function to initialize the internal data for the given configuration
    void init(const int * const config) {

      // populate the full pairing matrix
      for (int i = 0; i < _ns; i++)
        _pm[i*_ns+i] = pcps::zero<SCALAR>();
      for (int i = 0; i < _ns; i++)
      for (int j = i+1; j < _ns; j++) {
        const SCALAR x = _wfn.pairing_matrix()[pcps::cmpd_pair_index(i,j,_ns)];
        _pm[i*_ns+j] =  x;
        _pm[j*_ns+i] = -x;
      }

      // initialize the inverse of the occupied pairing matrix and get the absolute value of the initial pfaffian
      pcps::pfn::init_matrices(_np, _ns, config, _pm, _occ_sites, _row_indices, _inv_mat, _work_mat, _ipiv, _iwork);

      // store the wavefunction's value at this configuration
      _init_val = this->get_wfn_value(config);

      // set the initialized flag to true
      _initialized = true;

      // reset the update counter
      _update_count = 0;

    }

    ///////////////////////////////////////////////////////////////////////////////////////////////////
    /// \brief   computes the wavefunction value for the given configuration
    ///
    /// \param[in]      config     The lattice configuration whose wavefunction value is desired.
    ///
    /// \return the value of the wavefunction in mantissa/exponent format.
    ///
    ///////////////////////////////////////////////////////////////////////////////////////////////////
    pcps::MantExp<SCALAR> get_wfn_value(const int * const config) {

      // initialize the wavefunction value
      pcps::MantExp<SCALAR> retval(pcps::unity<SCALAR>());

      // get the correlators' contribution
      pcps::correlator_weight(_configs_per_site, _wfn.correlators(), config, retval);

      // get the pfaffian's contribution
      retval *= pcps::pfn::get_pfaffian_value(_np, _ns, config, _pm, _ipiv, _iwork, _work_mat);

      // return the wavefunction value
      return retval;

    }

    // function to return the value of the wavefunction coefficient for the configuration that the object was last initialized for
    SCALAR initial_value() const { return _init_val.get(); }

    // function for copying the state from another CoeffValueFunc
    void copy(const pcps::CoeffValueFunc<SCALAR, pcps::Part_Fermion, pcps::Ref_Pfaffian> & other) {
      if (this == &other) return;
      _initialized = other._initialized;
      std::memcpy( _pm         , other._pm         , sizeof(SCALAR) * size_t(_ns * _ns ) );
      std::memcpy( _inv_mat    , other._inv_mat    , sizeof(SCALAR) * size_t(_np * _np ) );
      std::memcpy( _work_mat   , other._work_mat   , sizeof(SCALAR) * size_t(_np * _np ) );
      std::memcpy( _wv1        , other._wv1        , sizeof(SCALAR) * size_t(_np       ) );
      std::memcpy( _wv2        , other._wv2        , sizeof(SCALAR) * size_t(_np       ) );
      std::memcpy( _wv3        , other._wv3        , sizeof(SCALAR) * size_t(_np       ) );
      std::memcpy( _row_indices, other._row_indices, sizeof(int   ) * size_t(_ns       ) );
      std::memcpy( _occ_sites  , other._occ_sites  , sizeof(int   ) * size_t(_np       ) );
      std::memcpy( _ipiv       , other._ipiv       , sizeof(int   ) * size_t(_np       ) );
      std::memcpy( _iwork      , other._iwork      , sizeof(int   ) * size_t(_ns       ) );
    }

    // function to get the coefficient ratio for two configurations (new_config in numberator), which differ at the sites in diff_sites
    SCALAR get_coeff_ratio(const int n_diff_sites,
                           const int * const diff_sites,
                           const int * const new_config,
                           const int * const old_config) const {

      // check that this object is initialized
      assert(_initialized);

      // check that the number of sites affected by the move is supported
      assert( n_diff_sites == 0 || n_diff_sites == 2 );

      // return 1 if no sites were changed
      if (n_diff_sites == 0)
        return pcps::unity<SCALAR>();

      // initialize the coefficient
      SCALAR retval = pcps::unity<SCALAR>();

      // compute the correlators' contribution to the new coefficient
      if (n_diff_sites == 2) {
        assert( diff_sites[0] % 2 == diff_sites[1] % 2 );
        retval *= pcps::correlator_weight(_configs_per_site, _wfn.touching_correlators(diff_sites[0], diff_sites[1]), new_config);
        retval /= pcps::correlator_weight(_configs_per_site, _wfn.touching_correlators(diff_sites[0], diff_sites[1]), old_config);
      //} else if (n_diff_sites == 4) {
      //  const int * const ds = diff_sites;
      //  retval *= pcps::correlator_weight(_configs_per_site, _wfn.touching_correlators(ds[0], ds[1], ds[2], ds[3]), new_config);
      //  retval /= pcps::correlator_weight(_configs_per_site, _wfn.touching_correlators(ds[0], ds[1], ds[2], ds[3]), old_config);
      } else throw pcps::Exception( (boost::format("%i is not a valid value for n_diff_sites in pcps::CoeffValueFunc::get_coeff_ratio") % n_diff_sites).str() );

      // multiply by the ratio of pfaffians
      retval *= pcps::pfn::pfn_ratio(_np, _ns, n_diff_sites, diff_sites, new_config, old_config, _pm, _inv_mat, _occ_sites, _row_indices, _wv1);

      // return the coefficient
      return retval;

    }

    // function to update the internal data for a new configuration, which differs from the old configuration at the sites in diff_sites
    void update(const int n_diff_sites,
                const int * const diff_sites,
                const int * const new_config,
                const int * const old_config) {

      // check that this object is initialized
      assert(_initialized);

      // check that the number of sites affected by the move is supported
      assert( n_diff_sites == 0 || n_diff_sites == 2 );

      // there is nothing to do if no sites were affected
      if (n_diff_sites == 0) return;

      // update the inverse agp matrix and the alpha and beta site lists
      pcps::pfn::update(_np, _ns, n_diff_sites, diff_sites, new_config, old_config, _pm, _inv_mat, _occ_sites, _row_indices, _wv1, _wv2, _wv3);

      // increment the update counter
      _update_count++;

      // if we have updated too many times without re-initializing, re-initialize
      if (_update_count > 1000)
        this->init(new_config);

    }

    // Function to get the ratio of the current coefficient's derivatives to the coefficient's value.
    // Derivatives with respect to each wavefunction variable are considered, although derivatives
    // with respect to most correlator elements are zero, so only non-trivial derivatives are returned.
    void get_compact_derivatives(const int * const config,
                                 int    * const cor_der_indices,
                                 SCALAR * const cor_der,
                                 SCALAR * const agp_der,
                                 SCALAR * const orb_der,
                                 SCALAR * const par_der) const {

      // compute correlator derivative ratios (skip fixed correlators)
      if (_get_cor_der) {
        int c = 0;
        for (typename std::vector<pcps::Correlator<SCALAR> >::const_iterator corr = _wfn.correlators().begin(); corr != _wfn.correlators().end(); corr++)
          if (!corr->fixed()) {
            int cmpd = 0;
            const int * sites = &corr->site_vec()[0];
            const int nsites = corr->site_vec().size();
            for (int i = 0; i < nsites; i++)
              cmpd = cmpd * _configs_per_site + config[sites[i]];
            cor_der_indices[c] = cmpd;
            cor_der[c] = pcps::unity<SCALAR>() / (*corr)[cmpd];
            c++;
          }
      }

      assert( !_get_agp_der );
      assert( !_get_orb_der );

      // compute derivative ratios for the pairing matrix
      if (_get_par_der)
        pcps::pfn::der_ratios(_np, _ns, _inv_mat, _occ_sites, _row_indices, par_der);

    }

  };

} // end namespace pcps

#endif
